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PREFACE 
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propulsion  test  facilities,  AEDC,  AFSC,  Arnold  Air  Force  Station,  TN  37389.  Dr.  Koenig 
contributed  to  the  effort  as  a  consultant  from  the  Department  of  Aerospace  Engineering, 
Mississippi  State  University,  Mississippi  State,  Mississippi.  The  manuscript  was  submitted 
for  publication  on  April  1,  1983. 
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1.0  INTRODUCTION 

The  governing  field  equations  which  describe  the  general  cases  of  fluid  flow  are  often 
referred  to  as  the  Navier-Stokes  equations.  If  the  fluid  is  assumed  to  be  inviscid,  then  the 
equations  of  motion  are  often  called  the  Euler  equations.  If  the  fluid  is  assumed  to  be 
inviscid  and  the  flow  is  assumed  to  be  irrotational,  then  the  flow  is  termed  potential.  The  full 
Navier-Stokes  equations  are  a  coupled  system  of  second  order,  nonlinear  partial  differential 
equations;  the  Euler  equations  are  a  coupled  system  of  first  order,  nonlinear  partial 
differential  equations,  and  the  governing  equation  for  potential  flow  is  a  second  order, 
nonlinear  partial  differential  equation.  The  Navier-Stokes  equations  can  treat  virtually  any 
flow  problem  but  are  the  most  difficult  to  solve.  The  potential  flow  equation,  although  the 
easiest  to  solve,  is  the  most  limited  in  its  applications.  The  Euler  equations  are  intermediate 
in  both  respects  and  provide  a  practical  means  of  solving  fairly  difficult  problems  in  an 
economical  way.  If  used  properly,  the  Euler  equations  may  be  used  to  describe  three- 
dimensional,  unsteady  flow  at  any  Mach  number  and  without  assumptions  of  irrotational  or 
isoenergetic  conditions.  The  Euler  equations  are  not  subject  to  irrotational  and  isoenergetic 
assumptions  with  the  important  consequence  that  flows  with  curved  shocks  or  multiple, 
intersecting  shocks  can  be  treated,  thus  making  the  equations  applicable  to  transonic  and 
supersonic  flow  past  complex  bodies  or  combinations  of  bodies. 

2.0  BASIC  CONCEPTS 

The  computer  program,  ARO-1  (Ref.  1),  solves  the  three-dimensional,  unsteady  Euler 
equations  in  Cartesian  coordinates  using  a  finite  volume  (volume  flux)  approach.  The  basic 
numerical  algorithm  is  the  explicit  predictor-corrector  scheme  of  MacCormack  (Ref.  2). 
Reference  1  discusses  in  detail  this  computer  program,  the  fluid  mechanics  of  the  problem, 
and  the  numerical  techniques  involved  in  the  solution  to  the  equations.  This  report  serves 
two  purposes:  (1)  to  introduce  the  principles  of  computational  fluid  dynamics  to  engineers 
and  scientists  unfamiliar  with  this  subject  and  (2)  to  provide  detailed  instructions  on  the  use 
of  this  particular  computer  solution  of  the  Euler  equations. 

2.1  THE  MESH 

The  Euler  equations  are  partial  differential  equations.  Numerical  solutions  to  differen¬ 
tial  equations  require  approximating  continuous  functions  (such  as  pressure)  and  con¬ 
tinuous  operations  (mainly  differentiation  and  integration)  by  discrete  values,  differences, 
and  sums  done  for  discrete  points  in  space  and  time.  The  discrete  points  in  space  form  the 
basis  of  what  is  called  a  computational  grid  or  mesh.  This  mesh  is  an  essential  part  of  any 
CFD  code.  It  must  be  fine  enough  to  resolve  important  flow  details,  yet  coarse  enough  to 
keep  storage  and  computational  time  requirements  reasonable.  In  addition,  it  should  be  laid 
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out  and  numbered  in  a  consistent  manner  so  that  the  computations  can  progress  smoothly 
and  unambiguously  through  the  mesh.  The  mesh  is  also  one  of  the  portions  of  this  program 
that  users  will  need  to  modify  or  create  on  their  own  for  particular  problems;  therefore,  a 
thorough  understanding  of  this  mesh  scheme  is  necessary. 

Basically  a  mesh  or  grid  should  concentrate  points  near  the  body  and  spread  them  out 
farther  away.  Figure  1  presents  an  example  of  a  very  coarse  grid  about  an  axisymmetric 
body.  Figure  la  shows  the  side  view,  Fig.  lb  shows  the  front  view  midway  back  of  the  body, 
and  Fig.  lc  shows  an  oblique  view  of  a  portion  of  the  grid.  Nomenclature  used  in  the 
program  is  also  given.  The  mesh  shown  here  would  be  referred  to  as  a  10  x  3  x  7  mesh;  i.e., 
there  are  10  planes  perpendicular  to  the  X  axis,  3  surfaces  enclosing  the  X  axis,  and  7  planes 
radially  out  from  the  X  axis.  The  numbering  scheme  increases  the  indices  with  increasing 
values  of  X,  radial  location,  and  azimuthal  location  (clockwise  from  vertical).  The 
subscripts  (JX,  JY,  JZ)  are  not  synonymous  with  the  coordinates  (X,  Y,  Z).  For  example,  in 
this  grid  two  points  with  the  same  JZ  may  have  different  Z  values. 

The  network  of  lines  forms  six-sided  volumes  (some,  on  the  X  axis,  degenerate  to  only 
five  sides)  with  eight  corners  or  vertices.  To  refer  to  a  particular  volume  in  the  mesh,  the 
subscripts  for  the  corner  in  the  most  downstream,  radially  outward,  clockwise  location  are 
used.  The  subscripts  for  this  corner  are  also  the  largest  values  of  JX,  JY,  and  JZ  which  any 
corner  in  the  volume  has.  For  example,  the  shaded  volume  in  Fig.  lc  would  be  designated  as 
volume  (5,  3,  2).  The  physical  dimensions  of  the  cell  volumes  and  surface  areas  are  referred 
to  as  metrics.  Determining  these  quantities  is  not  as  straightforward  as  might  be  expected; 
this  will  be  discussed  in  more  detail  elsewhere  in  this  report. 

In  ARO-1  the  grid  is  generated  (or  input)  in  one  subroutine,  CONE,  whereas  the  metrics 
are  determined  in  another  subroutine,  YOLGRD.  Usually  the  grid  is  generated  in  an 
independent  mesh  program,  and  CONE  is  used  to  transmit  the  coordinates  (X,  Y,  Z)  from 
disk  storage  to  ARO-1.  The  subroutine  VOLGRD  is  quite  general  and  will  work  for  any  grid 
provided  the  numbering  scheme  described  above  is  used. 

2.2  THE  EQUATIONS 

The  Euler  equations  describe  the  flow  of  an  inviscid  fluid.  For  negligible  body  forces 
(mainly  gravity)  these  equations  can  be  written  in  rectangular  Cartesian  coordinates  as 

Continuity  +  V  *  =  0  (1) 
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Momentum 


Energy 


d_ 

dt 


a(gu) 

at 

a(ev) 

at 

a(gw) 

at 


+  V  .  (euq)  =  -■ 
+  V  •  (evq)  =  - 
+  V  ■>  (gwq)  =  — 


dp_ 

ax 

dp 

aY 


V  •  (gq) 


(2) 

(3) 


where  q  =  ui  +  vj  +  wk.  In  addition,  an  equation  of  state  is  required.  Here  a  perfect  gas  is 
assumed,  so  that 

p  =  gRT 

e  =  cvT  =  internal  energy 

7  =  cp/cy 

c2  =  7RT  =  Yp/g  =  speed  of  sound  (squared)  (4) 


With  some  algebra  these  equations  can  be  rearranged  so  all  have  the  same  strong  con¬ 
servation  form: 


+  V  .  F  =  0  (5) 

at 

Here  G  is  a  vector  formed  of  the  dependent  flow  variables  g,  gu,  gv,  gw,  and  e;  and  F  is  a 
matrix  of  flux  quantities,  such  as  guv  (the  flux  of  X momentum  in  the  Y  direction).  More 
exactly,  it  is  expressed  as 


“e  " 

gu 

gv 

Q  W 

gu 

gu2  +  p 

guv 

Q  UW 

gv 

F  = 

guv 

gv2  +  p 

QV  W 

gw 

guw 

gvw 

Q\V2  +  p 

e 

_(e  +  p)u 

(e  +  p)v 

(e  +  P)wea 
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The  energy  e  here  is  actually  the  sum  of  the  internal  and  kinetic  energy  per  unit  volume: 


e  =  g^e  4 — q2  =  u2  4-  v2  4- 


(7) 


and  is  related  to  pressure  by 


P  1  ? 

e  =  +  T  eq2 


(8) 


If  Eq.  (5)  is  applied  to  one  of  the  small  stationary  volumes,  V,  of  the  mesh  it  may  be 
written  as  (see  Section  4.0) 


dG 

dt 


--LJ 

vs 


n  dS 


(9) 


where  n  is  the  outward  unit  normal  to  the  surface,  S,  enclosing  V,  and  G  is  the  average  value 
of  G  in  the  volume.  Here  this  average  value  is  approximated  by  the  value  at  the  center. 

Equation  (9)  is  the  (schematic)  governing  equation  of  Euler  flow.  It  describes  three- 
dimensional  unsteady  flows,  but  is  used  in  ARO-1  to  calculate  steady  flows.  To  do  so,  a  set 
of  initial  conditions  for  G  (and  therefore  F)  is  assumed.  These  initial  conditions  are  arbitrary 
but  imposing  free-stream  conditions  everywhere  works  well.  Through  the  use  of  Eq.  (9),  the 
time  rate  of  change  of  G  is  determined,  and  values  of  G  at  later  times  are  found.  Time  is 
allowed  to  continue  until  the  values  of  G  at  one  time  are  sufficiently  close  to  the  values  at  the 
preceding  time  to  consider  the  process  as  converged.  These  converged  values  of  G  are  then 
the  solution  to  the  problem.  The  solution  is,  of  course,  constrained  by  boundary  conditions. 
These  boundary  conditions  ensure  the  uniqueness  of  the  solution. 

The  ARO-1  program  solves  Eq.  (9)  using  MacCormack’s  explicit  predictor-corrector 
scheme.  Details  of  this  scheme  as  applied  to  Eq.  (9)  are  described  in  the  next  section.  Here  it 
is  sufficient  to  observe  a  few  general  characteristics  of  the  technique  and  this  equation. 
Equation  (9)  provides  a  means  of  finding  G  at  a  later  time  knowing  F  at  the  present  time. 
The  explicit  nature  of  the  technique  means  that  at  a  given  instant  of  time  all  values  of  F  are 
known  and  3G/3t  can  be  found  by  one  equation.  (Implicit  techniques  require  solving 
simultaneous  equations.)  Because  calculations  are  actually  done  using  finite  time  steps  and 
values  of  F  known  only  at  certain  points  (rather  than  continuous  distributions),  the  time 
derivative  and  surface  integrals  in  Eq.  (9)  can  only  be  approximated.  The  predictor-corrector 
nature  of  the  technique  means  that  Eq.  (9)  is  used  twice,  first  to  predict  values  of  G  and  F  at 
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at  a  later  time,  then  to  take  these  predicted  values  and  obtain  improved  or  corrected  values 
of  G  and  F  at  this  later  time.  These  corrected  values  can  be  shown  to  be  formally  second 
order  accurate  in  the  time  and  space  intervals,  which  is  sufficient  for  many  applications. 

3.0  IMPLEMENTATION  OF  MACCORMACK’S  EXPLICIT 
PREDICTOR-CORRECTOR  TECHNIQUE 

The  equations  of  motion  are  schematically  written 

J|L  +  V  »  F  =  0  (10) 

1 

where  G  represents  the  flow  variables  and  F  represents  the  flux  vectors.  If  G  were  known  at 
all  times  throughout  space,  the  problem  would  be  solved.  The  flow  variable  G  is  found  as 
follows. 

Suppose  at  some  time  t  we  know  G  everywhere  (either  because  we  assumed  some  values 
or  have  calculated  them  up  to  that  time).  The  values  of  G  at  a  later  instant  in  time  can  be 
found  by  expanding  G  in  a  Taylor  series  about  t,  i.e, 

G(t  +  At)  =  G(t)  +  At  +  ilV).  at2  +  •  •  •  <H) 

where  the  subscript  t  on  the  derivatives  means  to  evaluate  the  derivatives  at  time  t.  The  larger 
the  time  step  At  is,  the  more  terms  are  needed  in  Eq.  (11)  to  give  a  good  value  for  G(t  +  At). 
For  finite  time  steps,  Eq.  (11)  might  be  approximated  by 

G(t  +  At)  =  G(t)  +  (-£-)a>g  At  (12) 

where  (3G/dt)avg  is  some  effective,  or  average,  value  of  the  time  derivative  over  the  interval  t 
to  t  +  At.  Finding  the  new  value  of  G  now  becomes  a  problem  of  finding  (dG/dt)avg.  This  is 
done  through  the  use  of  Eq.  (10).  The  flux  vectors  F  in  Eq.  (10)  are  simply  algebraic 
combinations  of  the  flow  variables  G.  If  G  is  known  everywhere  at  time  t,  then  F  is  also 
known  everywhere  at  the  same  instant,  and  Eq.  (10)  gives  a  way  to  find  the  time  derivative  of 
G. 

To  use  Eq.  (10)  on  a  grid  of  discrete  points,  some  manipulations  are  necessary.  Equation 
(10)  is  rearranged  and  integrated  over  one  of  the  small,  stationary  grid  volumes  to  give 

I  ~~  dv  =  -  j  V  .  F  dV  (13) 

v  v 
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Consider  the  left-hand  side  of  Eq.  (13).  Since  the  volumes  are  assumed  to  be  stationary,  the 
time  derivative  can  be  taken  outside.  The  mean  value  theorem  of  calculus  can  then  be 
applied  to  the  integral.  This  sequence  yields 

l  irdV  -  i  l Gdv  -  £<°  l dv>  -  l<5v)  = <14> 

where  G  is  the  mean  value  of  G  in  the  volume  V.  If  the  volume  is  small,  then  G  is 
approximately  the  same  as  the  value  of  G  at  the  volume  center;  therefore,  Eq.  (14)  gives 

(  (jy  =  yJG_  (15) 

J  at  at 

V 

Note  that  Eqs.  (13)  and  (14)  are  exact  mathematically,  given  that  Eq.  (10)  and  Eq.  (11)  are 
also  exact  (with  the  proper  conditions  of  G  being  differentiable  and  integrable).  Equations 
(12)  and  (15)  introduce  approximations.  By  considering  the  right-hand  side  of  Eq.  (13),  we 
use  the  divergence  theorem  to  change  the  volume  integral  to  a  surface  integral,  i.e., 

j  V  •  F  dV  =  }  F  •  n  dS  (16) 

v  s 


where  h  is  the  unit  outward  normal  of  S.  We  use  Eqs.  (15)  and  (16)  in  Eq.  (13)  and  divide 
through  by  V  to  give 

-^--^iF.&dS  (17) 

s 

Equation  (17)  is  now  the  means  by  which  the  first  time  derivative  of  G,  evaluated  at  the 
volume  center  and  at  time  t,  can  be  found  approximately.  The  quantity  G  on  the  right-hand 
side  of  Eq.  (15)  and  in  Eq.  (17)  is  the  same  as  the  G  in  Section  2.2. 

To  eliminate  possible  confusion,  the  basic  purpose  must  be  reexamined.  Our  goal  is  to 
find  G(t  +  At)  knowing  G(t).  We  decided  to  do  this  by  Eq.  (12),  but  this  equation  requires 
finding  (3G/dt)avg  (which  has  has  not  been  defined  yet).  A  specific  time  derivative  can  be 
found  using  Eq.  (17);  therefore,  we  have  a  tool  to  use.  However,  there  is  an  additional 
consideration  to  using  this  tool.  Equation  (17)  requires  an  integration  of  F  over  the  surfaces 
of  a  grid  volume,  but  we  only  know  F  at  discrete  points  in  space,  in  fact,  only  at  the  centers 
of  volumes.  So  we  have  two  problems  to  solve  in  order  to  find  G(t  +  At):  (1)  how  do  we 
find  (3G/3t)avg  and  (2)  how  do  we  perform  the  integration  in  Eq.  (17)?  These  two  problems 
are  addressed  simultaneously  by  MacCormack’s  explicit  predictor-corrector  technique. 
Justifying  that  the  steps  in  MacCormack’s  method  are  valid  and  accurate  to  a  particular 
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order  is  not  that  simple,  but  it  can  be  done.  Let  us  accept  that  the  technique  does  work  and  is 
essentially  second  order  accurate  and  instead  worry  only  with  how  it  is  used  and  that  it,  at 
least,  is  reasonable. 

A  digression  is  necesary  at  this  point  to  review  the  computational  grid.  The  grid  points 
form  the  corners  of  six-sided  volumes  (hexahedra),  three  of  which  are  sketched  below. 


The  crosses  in  the  sketch  are  the  volume  centers  where  the  flow  variables  and  flux  vectors 
are  known,  here  denoted  by  F_,  F,  and  F+ .  The  two  shaded  faces,  Sk,  and  Sk,  are  opposing 
faces  of  the  middle  volume  to  serve  as  examples. 

The  integration  in  Eq.  (17)  determines  the  net  flux  through  the  volume  in  question.  Since 
each  volume  has  three  pairs  of  opposing  sides,  the  integral  can  be  treated  as  the  sum  of  the 
integrals  over  each  of  the  six  surfaces,  i.e., 

3  3 

(  F*  n  dS  =  E  (  F  •  n  dS  +  E  (  F  •  n  dS 

k=i  .  k=l  _ 

s  Sir  sk  (18) 

which  is  exact.  However,  F  is  not  known  over  each  surface  but  only  at  the  volume  center 
points,  so  the  individual  integrals  must  be  approximated.  This  approximation  can  be  done  in 
a  variety  of  ways.  For  example,  the  flux  across  Sk  might  be  determined  using  an  average  of 
F~  and  F,  or  an  average  of  F~ ,  F  plus  other  neighboring  centers  (not  shown)  or  using  F_ ,  or 
F  alone.  The  MacCormack  method  essentially  uses  an  average  of  F_  and  F,  but  does  so 
indirectly  in  two  steps,  as  will  be  shown  later  in  this  report. 

The  term  predictor-corrector  means  a  solution  for  G  (t  +  At)  is  first  predicted,  then  this 
prediction  is  corrected.  The  quantity  (dG/dt)avg,  defined  by 
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(4),,,  -  4- [(4), +  (4)]  (19) 

incorporates  this  idea;  (3G/3t)t  predicts  a  value  and  (dG/dt)  corrects  the  result.  (Note  that 
the  overbar  here  indicates  a  predictor-corrector  quantity  and  is  a  different  usage  from  that 
of  Section  2.2.)  The  predicted  value  of  G  (t  +  At)  is  given  by 

G(t  +  At)  =  G(t)  +  At  (20) 

Equations  (17)  and  (18)  are  used  to  obtain  (3G/dt)t: 

(^■)l--v(t?,F-sf+tl,F“-s0  (21) 

Note  here  that  the  flux  integrals  are  approximated  by  computing  the  scalar  product  of  the 
flux  vector  “behind”  the  appropriate  face  with  the  area  vector  of  that  face.  The  quantities 
G  (t  +  At)  from  Eq.  (20)  can  now  be  used  to  find  predicted  values  for  the  flux  vectors, 
F(t  +  At).  These  new  flux  vectors  are  used  to  find  the  predicted  derivative,  which  serves  as  a 
corrector  to  G(t  +  At): 


(4),  “  - v(  J,  •  Sk+  +  J,  F  •  sr)  (22) 

Here  the  flux  integrals  are  approximated  by  using  the  flux  vectors  “ahead  of”  the 
appropriate  face.  Equation  (22)  is  the  corrector  step  and  completes  the  flux  averaging 
problem  mentioned  earlier. 

In  principle  all  the  ingredients  to  find  G(t  +  At)  are  now  available.  First,  using  the 
known  conditions  at  time  t,  we  find  (3G/3t)t  from  Eq.  (21).  By  substituting  this  result  into 
Eq.  (20)  to  find  G(t  +  At),  we  then  use  G(t  +  At)  to  find  F(t  +  At);  and  substituting 
F(t  +  At)  into  Eq.  (22),  Eq.  (22)  into  Eq.  (19),  and  Eq.  (19)  into  Eq.  (12)  yields  G(t  +  At). 
In  the  notation  of  Ref.  1,  G(t)  =  Gn,  G(t  +  At)  =  Gn+1,  and  F(t  +  At)  =  F.  Except  for  a 
factor  w,  Eqs.  (20)  and  (21)  here  form  the  first  part  of  Eq.  (7)  in  Ref.  1,  whereas  Eqs.  (12), 
(19),  (20),  and  (22)  here  combine  to  form  the  second  part  of  Eq.  (7)  in  Ref.  1. 

4.0  DISCUSSION  AND  RECOMMENDATIONS 

Time-dependent,  predictor-corrector  solution  techniques  of  the  three-dimensional, 
unsteady  Euler  equations  can  be  effectively  applied  to  a  rather  wide  range  of  practical 
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situations,  as  demonstrated  in  Refs.  1,  3,  and  4.  The  preceding  sections  have  presented 
introductions  to  three  essential  aspects  of  this,  or  any,  CFD  program:  the  equations  being 
solved,  the  computational  grid  on  which  the  solution  is  determined,  and  the  numerical 
algorithm  used  to  solve  the  equations.  Satisfactory  application  of  this  program  requires 
users,  especially  those  new  to  CFD,  to  understand  thoroughly  the  details  of  the  code.  Also 
included  in  this  report  are  a  discussion  of  program  organization  (Appendix  A),  a  discussion 
of  program  sequence  (Appendix  B),  and  a  detailed  description  of  the  major  subroutines 
(Appendix  C).  Although  some  of  this  information  is  not  absolutely  necessary  in  order  to  use 
the  program,  careful  study  of  the  Appendices  will  enable  users  to  apply  the  code  more 
efficiently  to  their  particular  problem.  The  full  capabilities  of  this  program  can  be 
determined  and  exploited  only  if  the  code  is  tested  and  used  for  a  variety  of  problems; 
therefore,  careful  study  before  application  is  highly  recommended. 
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APPENDIX  A 

PROGRAM  ORGANIZATION 


1.  Introduction 

The  program  in  its  present  form  consists  of  a  main  calling  program  and  17  subroutines. 
The  calculations  involve  considerable  looping  among  the  main  program  and  the  other  17 
subroutines.  Extensive  use  is  made  of  labelled  common  regions  to  transfer  data  among  these 
subroutines.  Input  to  the  program  is  quite  minimal.  Only  four  data  cards  are  read  although 
this  is  a  little  misleading,  since  the  mesh  coordinates  are  usually  input  through  an 
unformatted  READ  statement  in  one  of  the  subroutines  and  other  information  is  provided 
through  the  use  of  PARAMETER  and  DATA  statements. 

2.  Subroutines 

A  list  of  the  main  program  and  subroutines  follows  with  very  brief  descriptions  of  the 
main  activities  of  each  unit.  The  list  is  in  the  order  in  which  the  routines  appear  in  the 
program  and  not  in  the  order  in  which  they  are  called.  Some  of  these  subroutines  will  be 
described  in  greater  detail  in  later  sections. 

a.  AROl  is  the  main  calling  program.  It  calls  (in  chronological  order)  the  sub¬ 
routines  which  (1)  allow  for  restart,  (2)  establish  the  initial  conditions,  (3) 
generate  (or  read)  the  grid,  (4)  compute  the  grid  metrics,  (5)  initialize  certain 
arrays,  (6)  print  the  grid  mesh  point  locations,  (7)  print  the  initial  conditions, 

(8)  control  the  computational  cycles,  and  (9)  stop  the  program  when  everything 
is  finished. 

b.  START  subroutine  reads  a  restart  file  produced  from  a  previous  calculation 
and  provides  for  creation  of  a  subsequent  restart  file  at  the  conclusion  of 
additional  calculations.  This  subroutine  also  initializes  the  flux  accumulators, 
the  boundary  conditions,  and  time  step. 

c.  PVAR  subroutine  is  user-adaptable  and  produces  the  desired  printed  output  of 
flow-field  and  geometry  variables.  The  user  can  tailor  the  subroutine  as 
required. 

d.  PARRAY  subroutine  writes  any  array  passed  to  it. 

e.  CONE  subroutine  is  the  vehicle  for  inputting  the  mesh  coordinates.  Usually  the 
coordinates  are  computed  using  a  separate  mesh  generation  code;  however,  for 
simple  geometries,  the  coordinates  could  be  computed  in  this  subroutine. 
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f.  ICEES  subroutine  computes  the  initial  conditions  throughout  the  flow  field, 
based  on  the  specified  Mach  number  and  attitude. 

g.  VOLGRD  subroutine  computes  the  grid  metrics  (projected  areas  and  volumes). 

h.  BC  subroutine  enforces  the  boundary  conditions  [Eq.  (9),  Ref.  1]  and  calls 
subroutine  INFLOW. 

i.  INFLOW  subroutine  enforces  the  inflow  boundary  conditions  using  a  reference 
plane  characteristic  approach  (see,  for  example,  Kneile  et  al.*)„  For  subsonic 
flow  the  subroutine  uses  the  specified  total  temperature,  total  pressure,  and 
flow  angles  to  compute  the  required  dependent  variables.  For  supersonic  flow, 
all  variables  are  specified  at  inflow.  INFLOW  calls  subroutine  CENTER. 

j.  CENTER  subroutine  takes  the  mesh  coordinates  (X,  Y,  Z),  which  are 
coordinates  of  the  volume  corners,  and  computes  coordinates  of  the  volume 
centers. 

k.  FLUXX,  FLUXY,  FLUXZ  are  subroutines  which  compute  the  fluxes  of  mass, 
momentum,  and  energy  through  the  faces  of  the  grid  volumes  in  the  JX,  JY, 
and  JZ  directions,  respectively.  Flux  accumulators  which  sum  the  net  flow  for 
the  predictor-corrector  calculations  are  also  computed.  Second  order  smooth¬ 
ing  of  the  dependent  variables  is  also  available,  depending  on  the  value  of  the 
parameter  KSMOOTH. 

l.  UPDATE  subroutine  computes  the  new  or  updated  values  of  the  flow  variables 
(G)  using  the  flux  accumulators.  This  is  where  parts  of  Eq.  (7)  of  Ref.  1  appear 
explicitly. 

m.  DTCAL  subroutine  computes  the  new  time  step  to  be  used  [Eq.  (8),  Ref.  1]. 

n.  STEP  subroutine  controls  the  calls  to  the  FLUX,  UPDATE,  SMOOTH,  and 
DTCAL  subroutines. 

o.  BSTEP  subroutine  controls  the  calls  to  the  STEP  and  PVAR  subroutines. 

p.  SMOOTH  subroutine  computes  smoothed  or  averaged  values  of  the  flow 
variables  inside  the  grid. 


*Kneile,  K.  R,,  Todd,  D.  C.,  and  Jacocks,  J.  L.  “Characteristic  Boundary  Conditions  for 
ARO-1.”  AEDC-TR-82-28 ,  May  1983. 
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3.  Inputs 

Information  is  transferred  into  the  program  in  three  ways:  through  READ,  DATA,  and 
PARAMETER  statements.  For  normal  operation  only  four  data  cards  are  read.  However,  if 
the  program  is  to  be  run  with  a  restart  capability,  then  additional  information  may  be  read 
from  a  specified  file. 

The  data  cards  will  be  discussed  in  the  order  in  which  they  are  encountered,  with 
definitions  and  typical  values  of  the  variables  being  read.  The  following  applies  to  cal¬ 
culations  without  a  restart. 

READ  (5,100)  KIN,  KOUT,  KSMOOTH,  TFLAG 

100  FORMAT  (313,  E24.0) 

This  read  is  in  subroutine  START. 

KIN  =  logical  unit  number  of  restart  file  to  be  =0 

read  (0  =  no  restart) 

KOUT  =  logical  unit  number  of  restart  file  to  be  =  0 

written  after  calculations  are  complete 

KSMOOTH  =  flag  to  turn  on  second  order  =  0 

smoothing  in  FLUX  subroutines  (0  = 
no  smoothing) 

TFLAG  =  maximum  time  value  permitted  if  in-  =  7.00 

ternal  clock  of  computer  is  used  to 
time  the  computation,  min 

READ  (5,100)  FACT,  FSMACH,  ALPHA,  BETA,  PHI 

100  FORMAT  (8E10.0) 

This  read  is  in  subroutine  ICEES. 

FACT  =  parameter  for  computational  =0.8 

stability  (CFL  number) 

FSMACH  =  free-stream  Mach  number  =  0.8500 
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ALPHA  =  free-stream  flow  angle,  deg 
BETA  =  free-stream  flow  angle,  deg 


PHI  =  free-stream  flow  angle,  deg 


=  0 

=  0  ...  no 
incidence 

=  -10.00  .  . 
10-deg 
angle  of 
attack 

=  0 


More  information  on  the  free-stream  flow  angles  is  available  in  Appendix  C,  ICEES  and 
Free-Stream  Conditions. 


READ  (5,100)  NB,  NT,  NV,  FACTX 
100  FORMAT  (315,  F20.0) 

This  read  is  in  the  main  routine,  AROl.  A  very  important  parameter  is  the  product 


NB*NT 

=  total  number  of  time  steps  to  be  use 

=  4000  (for 
example) 

NB 

=  number  of  printouts,  i.e.,  calls  to 
PVAR 

=  20 

NT 

=  (NB*NT)/NB,  must  be  a  multiple  of  8 

=  200 

NV 

=  number  of  times  the  boundary-layer 
subroutine  is  called  for  versions  of  the 

code  with  viscous/inviscid 
interaction  capability  (Ref.  3) 

=  0 

FACTX 

=  for  restart  runs,  a  new  value  for 
FACT;  otherwise  set  to  zero  (For 
restart  runs,  subroutine  ICEES  is  not 
called  and  the  value  of  FACT  is  taken 
from  the  restart  file;  thus,  to  change 
to  a  new  value,  FACTX  is  used.) 

=  0 

READ  (5,100) 

100  FORMAT  (2I3,E24.0) 

This  read  is  in  subroutine  START.  A  blank  card  must  be  the  fourth  card. 
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After  NB  calls  to  subroutine  BSTEP,  the  variables  NB,  NT,  NV,  FACTX  are  read  again.  If 
a  zero  value  for  NB  is  encountered,  execution  is  terminated.  If,  after  a  specified  number  of 
iterations,  it  is  desired  to  alter  the  Courant  number  or  the  amount  of  printout,  this  data  card 
can  be  used  for  that  purpose  by  indicating  the  new  values  of  NB,  NT,  NV,  and  FACTX 
desired,  where  FACTX  is  the  new  Courant  number. 

Besides  these  reads  from  TAPE5,  there  are  several  READ(KIN)  statements  for  restart 
which  appear  in  subroutine  START. 

A  second  form  of  input  is  the  DATA  statement.  The  occurrences  of  DATA  are  as 
follows: 

In  AROl: 

DATA  GAM,  GM1,  GS2,  RGM1  /.../,  where 


GAM  =  7,  GM1  =  7-1,  GS2  =  7/2,  RGM1  =  (7-I)-1. 

The  other  two  DATA  statements  set  initial  values  for  certain  counters  and  reference 
numbers. 

In  FLUXX,  FLUXY,  FLUXZ: 

DATA  KODE  /  .  .  ./  provides  a  way  to  cycle  around  the  eight  corners  of  a  grid 

volume.  These  should  not  be  changed. 

The  last  input  method  is  the  PARAMETER  statement.  It  sets  values  of  certain  quantities 
which  then  cannot  be  changed  by  any  later  program  statement.  It  must  appear  in  every 
subroutine  using  those  quantities.  Here  only  one  statement  is  used,  but  it  is  in  almost  all  of 
the  subroutines. 

PARAMETER  (NNX  =  80,  NNY  =  32,  NNZ  =  2,  .  .  .  ) 

This  particular  statement  is  for  dimensioning  arrays  for  an  80  x  32  x  2  grid.  Other  grid 
sizes  can  be  specified  by  replacing  these  numbers.  This  must  be  done  in  every  PARAMETER 
statement.  One  major  benefit  of  using  the  PARAMETER  statement  is  that  variables  may  be 
dimensioned  (in  DIMENSION  or  COMMON  statements,  for  example)  with  variable  sizes,  if 
these  sizes  are  first  specified  in  PARAMETER  statements.  This  saves  a  tremendous  amount 
of  work  in  this  program  when  grid  sizes  are  changed.  (The  PARAMETER  statement  is  a 
feature  of  Cray  computers  and  will  require  replacement  when  ARO-1  is  installed  on  a 
different  computer.) 


19 


AEDC-TR-83-22 


4.  Variables 

There  are  large  numbers  of  variables  in  this  program,  many  of  them  triply  subscripted 
with  dimensions  of  about  (NNX+1,  NNY+1,  NNZ+1).  This  leads  to  large  storage 
requirements.  Most  of  these  variables  appear  in  labelled  COMMON  statements,  and  these 
labelled  COMMON  statements  appear  in  most  of  the  subroutines  even  though  the 
subroutines  may  actually  require  only  a  few  of  the  variables.  The  variable  names  are,  in 
many  cases,  descriptive  of  the  actual  quantities  they  represent,  thus  helping  in  understanding 
and  using  the  program.  The  following  list  gives  the  more  important  variables,  their 
dimension  size,  definitions,  and  the  subroutines  in  which  they  are  actually  used.  This  list 
does  not  include  some  insignificant  variables  and  variables  which  appear  exclusively  in 
SMOOTH.  (In  the  following,  the  designation  NNX1  means  NNX  +  1,  with  analogous 
meanings  for  NNY1  and  NNZI.) 


A  (NNX1,  NNY1,  NNZI) 

-  dummy  name  for  any  of  the  flow 
variables 

PARRAY 

ADA 

=  A  •  A,  magnitude  of  area  vector 
squared 

BC,  DTCAL 

ALPHA 

=  free-stream  flow  angle 

ICEES 

AX, AY,  AZ 

=  average  value  of  the  x,  y,  and  z 
components  of  the  area  vectors 
of  opposing  faces 

VOLGRD 

AX1,  AY1,  AZ1 

(NNX,  NNY1,  NNZI) 

=  components  of  area  vector  of 
face  with  constant  JX 

VOLGRD, 
FLUXX, 
DTCAL, 
INFLOW,  BC 

AX2,  AY2,  AZ2 

(NNX1 ,  NNY,  NNZI) 

=  components  of  area  vector  of 
face  with  constant  JY 

VOLGRD, 
FLUXY, 
DTCAL,  BC 

AX3,  AY3,  AZ3 

(NNX1,  NNY1,  NNZ) 

=  components  of  area  vector  of 
face  with  constant  JZ 

VOLGRD, 
FLUXZ, 
DTCAL,  BC 
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BETA 

CLSTEP 

DE,  DR,  DRU,  DRV,  DRW 
(NNX1,  NNY1,  NNZ1) 

DT,  DTM 
DTSTEP 

DTSVOL 

DUM 

DX1,  DY1,  DZ1 
DX2,  DY2,  DZ2 

DXC,  DYC,  DZC 

E  (NNX1,  NNY1,  NNZ1) 

EINIT 

FACT,  FACTX 

FF 


free-stream  flow  angle,  negative 
of  angle  of  attack 

coefficient  of  flux  summation 
[Eq.  (11),  Ref.  1] 


flux  accumulators  for  energy, 
density,  and  (X,  Y,  Z) 
momentum 


time  step  [Eq.  (8),  Ref.  1] 

coefficient  of  flux  accumulators 
[Eq.  (11),  Ref.  1] 

coefficient  of  flux  accumulators 

denominator  of  Eq.  (11),  Ref.  1 

components  of  one  face  diagonal 

components  of  opposite  face 
diagonal 

components  of  vector  between 
centers  of  opposing  faces 

energy 


initial  value  of  energy 
CFL  number 


2*  (gq  •  A)/(A  •  A),  mirror 
condition  [Eq.  (9),  Ref.  1] 


ICEES 


UPDATE 


START, 

SMOOTH, 

FLUX, 

UPDATE 

DTCAL 

UPDATE 

DTCAL 

UPDATE 

DTCAL 

VOLGRD 

VOLGRD 


VOLGRD 


ICEES,  BC, 

FLUX, 

DTCAL 

ICEES 

AROl, 

ICEES, 

DTCAL 

BC 
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FSMACH 

— 

free-stream  Mach  number 

FI,  F2,  F3,  F4,  F5 
(NNX1,  NNY1) 

= 

flux  of  mass,  (X,  Y,  Z) 
momentum,  and  energy 

GAM 

= 

7,  ratio  of  specific  heats 

GM1 

= 

7  -  1 

GS2 

= 

7/2 

INC 

= 

counter,  number  of  times  STEP 
is  called 

JJ,  JS,  JT 

= 

indices  for  predictor-corrector 

JX,  JY,  JZ 

= 

primary  indices 

JX1,  JX2,  JY1,  JY2 

JZ1,  JZ2 

= 

indices  for  predictor-corrector 

KIN 

= 

logical  unit  number  for  restart 
input  file 

KODE 

= 

code  for  predictor-corrector 
sequence 

KOUNT 

= 

counter  for  time  steps 

KOUT 

= 

logical  unit  number  for  restart 
ouput  file 

LIMIT 

= 

counter  limit  for  calls  to  STEP 

NB 

= 

number  of  printouts 

NEX,  NEY,  NEZ 

== 

grid  size  - 1 

PVAR 

ICEES 

FLUXX, 

FLUXY, 

FLUXZ 

ICEES 

PVAR 

ICEES 

DTCAL 

STEP 


FLUX 

various  places 
FLUX 

AROl 

START 

FLUX 

PVAR,  BC, 
STEP 

AROl 

START 

STEP 

AROl 

various  places 
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NNX,  NNY,  NNZ 

—  grid  size 

various  places 

NNX1,  NNY1,  NNZ1 

=  grid  size  +  1 

various  places 

NT 

=  frequency  of  printouts 

AROl 

BSTEP 

NV 

=  number  of  calls  to  boundary- 
layer  subroutine 

AROl 

OMEGA 

=  artificial  viscosity  parameter 
of  Eq.  (7),  Ref.  1 

UPDATE 

OMTO 

=  coefficient  of  flux  accumulators 
(1  -  2  «) 

UPDATE 

P  (NNX1,  NNY1,  NNZ1) 

=  pressure 

PVAR,  BC, 

FLUX, 

START, 

UPDATE, 

DTCAL, 

SMOOTH 

PHI 

=  free-stream  flow  angle 

ICEES 

PI 

=  7 r 

CONE 

PINIT 

=  initial  value  of  pressure 

ICEES 

QDA 

=  q  •  A,  volume  flux 

BC 

FLUX 

R  (NNX1,  NNY1,  NNZ1) 

=  density 

START, 

ICEES, 

BC,  FLUX, 

UPDATE, 

DTCAL, 

SMOOTH 

RGM1 

=  (7  "  I)"1 

DTCAL 
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RQDA 
RQDQ 

RQDRQ  =  e  Q  *  6  Q 

RU,  RV,  RW 

(NNX1,  NNY1,  NNZ1)  =  (X,  Y,  Z)  momentum 


=  eq  •  A 

=  eq  •  q 


RUINIT,  RVINIT, 

RWINIT 

TFLAG 

TIME 

VOL  (NNX1,  NNY1,  NNZ1) 

VOLMIN 
X,  Y,  Z 

(NNX1,  NNY1,  NNZ1) 

XC,  YC,  ZC  (NNX) 

XDUM 


=  initial  (X,  Y,  Z)  momentum 
=  internal  program  time  limit 

=  current  flow  time 

=  grid  cell  volume 

=  minimum  grid  cell  volume 

=  grid  coordinates 

=  coordinates  of  face  centers 

=  denominator  of  Eq.  (11), 
Ref.  1 


DTCAL 

SMOOTH, 

UPDATE, 

DTCAL 

DTCAL 


START, 

ICEES, 

BC,  FLUX, 
UPDATE, 
DTCAL, 
SMOOTH 


ICEES 

START 

STEP 

PVAR, 

PARRAY, 

UPDATE 

VOLGRD, 

UPDATE, 

DTCAL, 

SMOOTH 

VOLGRD 


PVAR, 

CONE, 

VOLGRD 

VOLGRD 

DTCAL 
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The  major  calls,  loops,  and  computations  (in  the  order  they  occur)  are  described  below. 
This  sequence  is  for  a  run  without  the  boundary-layer  calculations,  the  internal  clock,  or  the 
restart  option,  and  is  for  the  analysis  of  one  body  started  impulsively  with  free-stream 
conditions  everywhere.  The  activities  described  are  performed  in  the  program  or  subroutine 
identified  immediately  above  the  description. 

The  input  cards  will  be  taken  as 


KIN 

= 

0 

no  restart 

KOUT 

= 

0 

TFLAG 

= 

7.00 

FACT 

= 

0.8 

FSMACH 

= 

0.85 

ALPHA 

= 

0.0 

BETA 

= 

0.0 

PHI 

= 

0.0 

NB 

= 

20  | 

4000  times  steps  with  a  printout  every 

NT 

= 

200  j 

200  time  steps 

NV 

= 

0 

no  call  to  boundary-layer  routine 

FACTX 

— 

0.0 

These  are  typical  values  for  a  run  without  the  boundary  layer  or  restart. 

AROl 

calls  START 

reads  KIN,  KOUT,  KSMOOTH,  TFLAG 
if  KIN  =  0,  returns  to  AROl 
calls  ICEES 

reads  FACT,  FSMACH,  ALPHA,  BETA,  PHI 
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computes  initial  values  of  R,  RU,  RV,  RW,  E 
returns  to  AROl 
calls  CONE 


grid  points  (X,  Y,  Z)  input 
returns  to  AROl 
calls  VOLGRD 

computes  AX1,  AY1,  AZ1,  AX2,  AY2,  AZ2,  AX3,  AY3,  AZ3,  YOL 

returns  to  AROl 

calls  entry  MIDDLE  in  START 

computes  P 

initializes  flux  accumulators  DR,  DRU,  DRV,  DRW,  DE 
calls  BC 

computes  initial  values  of  boundary  conditions 
returns  to  MIDDLE 
calls  DTCAL 
initializes  time  step 

if  KIN  =  0  and  KOUT  =  0,  returns  to  AROl 
calls  PVAR 

writes  initial  flow  variables 

returns  to  AROl 

reads  NB,  NT,  NY,  FACTX 

At  this  point  a  nested  loop  process  begins  which  (for  the  most  general  situation)  is 
controlled  by  the  values  of  NB,  NT,  NV,  and  TFLAG.  For  the  sequence  described  here 
where  the  internal  program  clock  and  boundary  layer  are  not  being  used,  only  NB  and  NT 
control  the  process. 
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AROl  calls  BSTEP  (NB  times) 

BSTEP  calls  STEP  (NT  times) 

STEP  increments  KOUNT  by  1 

STEP  calls  FLUX  and  UPDATE  (2  times  for  each) 

FLUX  computes  DR,  DRU,  DRV,  DRW,  DE 

UPDATE  computes  TIME,  R,  RU,  RV,  RW,  E,  P, 

DR,  DRU,  DRV,  DRW,  DE 

UPDATE  calls  BC 

BC  computes  R,  E,  P,  RU,  RV,  RW  at  boundaries 

STEP  calls  SMOOTH  after  every  16  calls  to  FLUX  and  UPDATE 

SMOOTH  computes  DR,  DRU,  DRV,  DRW,  DE, 

R,  RU,  RV,  RW,  E,  P 

SMOOTH  calls  BC 

BC  computes  R,  E,  P,  RU,  RV,  RW  at  boundaries 
STEP  calls  DTCAL  after  every  128  calls  to  FLUX  and  UPDATE 
DTCAL  computes  DTSTEP 
BSTEP  calls  PVAR 

PVAR  writes  computed  flow  variables 

When  this  loop  sequence  is  finished,  KOUNT  will  equal  the  product  NB*NT;  there  will 
be  NB  printouts,  not  including  the  initial  flow  field,  and  the  control  will  be  back  in  AROl . 
For  the  case  of  no  restart,  the  remaining  steps  end  the  program. 

AROl  calls  entry  FINISH  in  START 
if  KOUT  =  0,  returns  to  AROl 
STOP 
END 
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APPENDIX  C 
SUBROUTINES 


1.  CONE 

This  subroutine  reads  the  (X,  Y,  Z)  coordinates  of  a  grid.  The  grid  is  produced  in  a 
separate  mesh  generation  program.  The  mesh  generation  routines  can  be  included  in  CONE 
if  desired. 

2.  ICEES  and  Free-Stream  Conditions 

Reference  1  states  that  all  variables  are  dimensionless,  with  the  reference  conditions 
usually  taken  to  be  free-stream  density  and  sound  speed,  the  latter  given  by 

-nr 

The  FORTRAN  program  variables  are  R,  P,  E,  RU,  RV,  and  RW.  Since  these  are 
dimensionless,  then,  for  example, 


and,  therefore, 


Similarly,  if  C  is  the  nondimensional  sound  speed,  then 


and,  therefore, 


=  1 


Since  and  cM  are  the  reference  values  (by  definition),  then 


P  = 


P 

6ooc2» 
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and  not  p/p^  as  one  might  guess.  But  from  Eq.  (4),  Ref.  1,  p  =  1/7  ec2,  so 

P  =  —  —  £c2  =  —  RC2 

7  Q  7 


and 


l_ 

7 


From  Eq.  (3),  Ref.  1,  we  have 


e  =  — - - h  —  gq2  =  —2—  +  —  7PM2 

7-1  2  7-1  2 


E  = 


e 


eooc 


2 

oo 


J — e_  + 

'y  —  1  o  2  p 

/  *  *=■  oo  oo  ^oo^oo 


and 


E-  =  (A)  Y  +  1  Ml 


,7-1/  7 

Summarizing  these  results,  the  free-stream  nondimensional  flow  is  decribed  by 


R»  =  J.c.  =  i.p.  =  f.E.  =  (^f  +  T 


Ml 


If  we  let  q  be  the  dimensional  velocity  vector  and  q  be  its  magnitude,  then 


Also, 


q  =  (u2  +  v2  +  w2)l/2 


U  = 


V 


-,  w  = 


w 


so  that  Ua,  Va,  and  are  simply  the  (X,  Y,  Z)  components  of  M^.  This  completes  the 
description  of  the  free  stream. 
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The  initial  conditions  are  taken  to  be  free-stream  conditions  everywhere  and  are 
established  in  ICEES.  There  are  five  equations  of  interest  in  ICEES  which  prescribe  initial 
values  of  RU,  RV,  RW,  P,  and  E,  i.e., 

RUINIT  =  Rm  =  X  component  of  Mk  (since  R^  =  1) 

RVINIT  =  =  Y  component  of  MB 

RWINIT  =  =  Z  component  of 

EINIT  =  Eto  =  (l/7-l)  1/7  +  1/2 

PINIT  =  Pw  =  I/7 

The  components  of  are  found  from  an  Euler  angle-type  approach  (see,  for  example, 
Etkin*).  The  orientation  is  slightly  unconventional,  however.  The  coordinate  frame  remains 
fixed  to  the  body,  and  the  incoming  flow  is  put  at  an  angle.  The  angles  a,  /3,  <j>  refer  to  yaw, 
pitch,  and  roll,  respectively,  where  a  is  the  angle  in  the  XZ-plane,  f3  is  the  angle  in  the  XY- 
plane,  and  4>  is  the  angle  in  the  YZ-plane.  For  the  coordinate  frame  as  established  in  Fig.  1,  a 
positive  angle  of  attack  (nose  up)  is  given  by  BETA  equal  to  a  negative  angle.  All  angles  in 
ICEES  are  in  degrees. 

3.  YOLGRD 

This  subroutine  computes  the  area  vectors  (or  project  areas)  of  the  grid  volume  faces  and 
the  volume  of  the  elements.  Before  the  actual  details  of  the  calculations  are  described,  a 
brief  discussion  of  the  philosophy  of  the  scheme  is  in  order. 


This  subroutine  was  written  to  be  able  to  handle  almost  any  conceivable  grid,  as  long  as 
the  grid  is  consistently  numbered.  In  many  cases  (such  as  in  this  program),  the  grid  elements 
are  taken  from  some  sort  of  regular  geometric  shape  with  well-defined  flat  sides;  however, 
this  is  not  always  the  case.  It  is  quite  possible  to  form  a  grid  such  that  the  corners  of  a 
volume  side  are  not  all  in  the  same  plane,  and,  consequently,  the  surface  is  not  flat,  as 
depicted  in  Fig.  C-l.  This  introduces  an  essential  difficulty.  When  the  side  is  flat,  its  area 
and  normal  vector  (which  together  give  the  area  vector)  can  be  found  unambiguously  and 
exactly  (assuming  the  corners  are  connected  by  straight  lines).  But,  if  the  four  corners  of  a 


*Etkin,  B.  Dynamics  of  Atmospheric  Flight.  John  Wiley  &  Sons,  Inc.,  New  York,  1972,  pp.  112-117. 
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side  are  not  coplanar,  then  the  area  is  not  well-defined  since  only  the  four  corners  are 
actually  known.  Even  if  the  surface  can  be  analytically  described,  it  will  have  a  normal 
vector  which  varies  in  direction  over  the  surface.  For  a  numerical  solution,  it  is  highly 
desirable  to  have  each  surface  of  a  volume  element  defined  by  a  single,  constant  area  vector. 
This  means  the  nonplanar  surface  must  be  represented  by  an  effective  area  vector  which 
must  be  determined  solely  in  terms  of  the  corner  coordinates.  In  addition,  the  method  for 
doing  this  should  yield  the  exact  answer  when  the  surface  is  flat.  Finally,  the  method  should 
be  capable  of  properly  handling  surfaces  that  degenerate  to  a  line  or  point  (as  happens  in 
almost  any  grid  scheme).  The  scheme  used  in  VOLGRD  satisfies  all  these  requirements. 


Figure  C-l.  Volume  element  with  nonplanar  surfaces. 

Determining  the  volume  of  an  element  has  similar  difficulties.  The  volume  of  a  flat-sided 
element  can  be  found  unambiguously  and  exactly  (although  the  calculations  may  be  quite 
tedious).  For  nonplanar  surfaces,  the  element  is  poorly  defined,  and  the  volume  becomes 
difficult  to  compute.  For  numerical  solutions,  this  computation  needs  to  be  simplified  and 
based  only  on  the  corner  coordinates.  Also,  the  result  must  be  exact  for  flat-sided  elements. 
The  technique  in  VOLGRD  satisfies  these  requirements  as  well. 
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A  representative  volume  is  shown  in  Fig.  C-2.  With  the  corners  labelled  as  shown,  this 
volume  would  be  referred  to  by  the  subscripts  (JX,  JY,  JZ).  The  surfaces  of  the  element  are 
surfaces  on  which  the  values  of  JX,  JY,  or  JZ  are  constant  and  are  denoted  as  Al,  A2,  and 
A3  respectively,  there  being  two  of  each  of  these  for  any  element.  The  two  Al  surfaces  have 
been  singled  out  here  to  serve  as  examples.  The  calculations  to  be  described  are  repeated  for 
the  other  two  surface  pairs.  All  surfaces  are  denoted  by  the  corner  having  the  largest 
subscripts;  therefore,  in  this  example,  the  two  Al  surfaces  are  referred  to  by  (JX-1,  JY,  JZ) 
on  the  left  and  (JX,  JY,  JZ)  on  the  right. 

Each  surface  has  two  diagonals  which  can  be  treated  as  vectors .  For  surface  Al  (JX,  JY, 
JZ),  the  diagonal  vectors  are 

D1  =  [X(JX,JY,JZ)  -  X(JX,JY-l,JZ-l)]i  +  [Y(JX,JY,JZ) 

-  Y(JX,JY-l,JZ-l)]j  +  [Z(JX,JY,JZ)  -  Z(JX,JY-  1,JZ-  l)]k 
=  DXli  +  DYlj  +  DZlk 


D2  =  [X(JX,JY-  1,JZ)  -  X(JX,JY,JZ-  l)]i  +  [Y(JX,JY- 1,JZ) 

-  Y(JX,JY,JZ-l)]j‘  +  [Z(JX,JY—  1,JZ)  -  Z(JX,JY,JZ-  l)]k 
=  DX2i  +  DY2j~  +  DZ2k 


Figure  C-2.  Metrics  of  a  volume  element. 
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The  area  vector  of  surface  A1(JX,  JY,  JZ)  is  defined  to  be  one-half  the  cross  product  of  the 
two  diagonal  vectors,  i.e., 


A1(JX,JY,JZ)  =  !/2(D1xD2)  =  !/2(DYl*DZ2  -  DZl*DY2)i  +  (DZ1*DX2- DXl*DZ2)j 

+  (DX1  *DY2  -  DY 1  *DX2)k  =  AXli  +  AYljV  AZlk 

which  is  exact  for  flat  surfaces  with  straight  lines  connecting  the  corners.  This  vector  is  also 
shown  in  Fig.  C-2.  Similarly  the  diagonal  vectors  and  area  vector  for  side  A1(JX-1,  JY,  JZ) 
are  found. 

Each  side  also  has  a  center  which  is  defined  to  be  the  average  of  the  corner  coordinates. 
For  example,  the  center  of  A1(JX,  JY,  JZ)  has  an  X  value  of 

XC(JX,JY,JZ)  =  !4[X(JX,JY,JZ)  +  X(JX,JY,JZ- 1)  +  X(JX,JY- 1,JZ- 1) 

+  X(JX,JY-  1 ,  JZ)] 

The  centers  of  opposing  surfaces  (here  A1(JX-1,  JY,  JZ)  and  A1(JX,  JY,  JZ))  can  be 
connected  by  a  line  which  also  can  be  expressed  as  a  vector: 

\  / 

,JY,JZ)  -  XC(JX— 1,JY,JZ)]?  +  [YC(JX,JY,JZ)  -  1jc(jx- 1,JY, 

+  [ZC(JX,JY,JZ)  -  ZC(JX-  l,JY,JZ)]k 

=  DXa  +  DYCj  +  DZCk 

The  two  area  vectors  of  opposing  faces  and  the  vector  between  the  centers  of  these  faces 
can  be  used  to  compute  a  volume  (not  necessarily  the  true  volume  of  the  element).  For  the 
example  given  here,  this  calculation  is 

VOL  =  Vi [A1(JX, JY, JZ)  +  A1(JX-  1,JY,JZ)]  •  DC1 

i.e.,  the  dot  product  between  the  average  of  the  opposing  face  area  vectors  and  the  vector 
connecting  the  centers  of  these  faces  gives  a  volume.  For  a  rectangular  parallelepiped,  this 
calculation  would  indeed  be  the  true  volume.  Any  pair  of  opposing  faces  could  be  used  for  a 
rectangular  parallelepiped,  and  the  result  would  be  the  same;  or  the  average  of  this  process 
for  the  three  pairs  could  be  used  to  recover  the  true  volume  of  a  rectangular  parallelepiped. 
For  an  arbitrary  element,  the  volume  is  defined  by  the  average  for  the  three  pairs,  ie., 


DC1  =  [XC(JX 


JZ)]j 
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VOL  -  y3{  [Al(JX)  +  Al(JX-l)]  «  DC1  +  [A2(JY)  +  A2(JY-1)]  •  DC2 
+  [A3(JZ)  +  A3(JZ-  1)]  •  DC3) 

where  the  subscript  notations  for  the  area  vectors  have  been  shortened. 

These  definitions  of  area  vectors  and  volumes  are  certainly  not  exact  for  arbitrary 
elements,  but  are  not  too  far  off  for  small  volumes.  For  regular  elements,  they  are  exact  and 
they  are  computationally  efficient.  These  schemes  are  similar  to  those  used  in  many  modern 
numerical  computational  programs. 

A  final  comment  on  the  nomenclature  used  in  the  preceding  discussion  is  necessary.  The 
equations  appearing  in  VOLGRD  are  in  terms  of  the  components  of  the  various  vectors;  the 
vectors  themselves  never  occur.  For  example,  the  symbol  D1  has  been  used  in  this  discussion 
for  one  of  the  diagonals.  In  VOLGRD  only  the  components  of  Dl,  denoted  DX1,  DY1, 
DZ1,  are  used.  Similarly,  A1  never  appears  but  AX1,  AY1,  and  AZ1  are  used,  and  so  on. 
Hopefully,  this  difference  in  notation  will  not  cause  any  confusion. 

4.  BC 

This  subroutine  enforces  the  boundary  conditions  on  the  computed  flow.  Typical 
boundary  conditions  for  a  body  in  an  unbounded  free  stream  are  free-stream  conditions  far 
away  from  the  body  and  a  specified  normal  velocity  at  the  body  surface  (usually  zero). 
Implementing  conditions  like  these  on  a  computational  mesh  of  finite  extent  requires  that 
these  conditions  be  handled  in  a  special  way.  The  technique  used  here  is  based  on  the 
concept  of  phantom  points,  i.e.,  points  which  are  physically  outside  of  the  computational 
mesh. 

Phantom  points  are  the  centers  of  volume  elements  which  lie  just  outside  of  the 
computational  mesh.  Figure  1  shows  a  (coarse)  grid  for  an  axisymmetric  body  which  would 
be  useful  for  zero  incidence  or  angle-of-attack  problems  (hence,  the  JZ  lines  go  halfway 
around  the  body).  The  first  volume  elements  in  the  mesh  to  be  encountered  as  one  proceeds 
downstream  are  those  having  JX  =  2,  since  the  volume  element  is  designated  by  a 
downstream  corner  (see  Fig.  C-3a).  If  one  referred  to  a  volume  (1,  JY,  JZ)  this  would  be 
immediately  upstream  of  the  mesh.  The  center  of  this  volume  is  a  phantom  point.  Now, 
referring  to  Figs.  C-3b  and  C-3e,  a  volume  denoted  by  (JX,  JY,  1)  would  be  just  to  the  left 
(or  counter  clockwise)  of  the  grid  volumes  above  the  body,  since  the  volume  element  is 
designated  by  a  clockwise  corner.  The  center  of  this  volume  is  also  a  phantom  point.  In  a 
similar  manner,  the  centers  of  volumes  (JX,  JY,  8),  where  NNZ  +1  =  8,  are  phantom 
points  just  clockwise  of  the  lower  grid  volumes  (Figs.  C-3b  and  C-3d.).  The  centers  of 
volumes  (NNX+1,  JY,  JZ)  are  phantom  points  just  downstream  of  the  last  streamwise 
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volumes  (Fig.  la).  Finally,  the  centers  of  volumes  (JX,  1 ,  JZ)  are  phantom  points  just  inside 
the  body,  since  volumes  are  designated  by  radially  outward  corners  (Fig.  C-3c). 


Figure  C-3.  Volume  elements  and  phantom  points. 
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The  phantom  points  are  used  in  the  following  way.  At  planes  of  symmetry  (for  example, 
lines  JZ  =  1  and  JZ  =  7  in  Fig.  1)  and  body  surfaces  (JY  =  1),  mirror  conditions  are 
applied.  Normally  on  a  solid  body  surface  or  on  a  plane  of  symmetry,  we  would  simply 
require  that  q  •  n  =  0  (q  =  velocity  vector,  n  =  normal  to  surface).  However,  in  the 
formulation  of  this  program  we  do  not  know  conditions  on  surfaces  but  only  at  the  centers 
of  volume  elements.  Consequently,  we  only  know  conditions  slightly  off  the  body  or  to  one 
side  of  a  plane  of  symmetry.  Mirror  conditions  dictate  that  the  velocity  vector  at  the 
phantom  point  inside  the  body  surface  or  on  the  other  side  of  the  symmetry  plane  has  a 
normal  component  equal  and  opposite  to  the  normal  component  of  velocity  at  the 
neighboring  true  volume  center.  One  might  say  that  the  phantom  point  velocity  vectors  are 
“mirror  images”  of  the  neighboring  true  velocity  vectors,  hence  the  name  mirror  conditions. 
Then,  for  small  volume  elements  these  conditions  essentially  impose  a  zero  normal  velocity 
on  the  surface  in  question.  Implementing  the  mirror  condition  leads  to  Eq.  (9)  of  Ref.  1: 


Qphantom  point  —  qtrue 


2  (—  ^ue)s 

V  S  .  S  / 


where  S  is  the  area  vector  of  the  surface  in  question.  In  subroutine  BC,  these  mirror 
conditions  are  imposed  at  the  phantom  points.  The  density,  energy,  and  pressure  at  these 
same  phantom  points  are  determined  by  imposing  a  zero  normal  gradient  condition;  i.e.,  R, 
E,  and  P  do  not  change  from  the  phantom  point  to  the  neighboring  volume  center.  Thus,  we 
have  R(JX,  1,  JZ)  =  R(JX,  2,  JZ),  R(JZ,  JY,  1)  =  R(JX,  JY,  2),  and  R(JX,  JY,  NNZ+  1) 
=  R(JX,  JY,  NNZ),  and  so  on. 

At  the  downstream  boundary  (JZ  =  NNX)  the  zero  normal  gradient  condition  is 
imposed  on  all  the  flow  quantities  R,  E,  P,  and  RU,  RY,  RW.  Thus,  we  have  R(NNX+  1, 
JY,  JZ)  =  R(NNZ,  JY,  JZ)  and  so  on.  This  condition  may  not  be  physically  true,  especially 
if  the  downstream  boundary  is  not  far  away  from  the  body,  but  it  has  been  found  to  work 
virtually  as  well  as  more  correct  (and  much  more  complicated)  outflow  conditions  and, 
therefore,  has  been  retained  in  the  program. 

5.  UPDATE 

This  subroutine  essentially  evaluates  Eq.  (7)  of  Ref.  1  [more  appropriately,  its  alternate 
form,  Eq.  (11)].  This  equation,  in  two  parts,  is  the  net  result  of  applying  the  predictor- 
corrector  technique  to  this  problem.  In  a  very  schematic  way,  both  parts  of  Eq.  7  [or  Eq. 
(11)]  look  like 

Gnew  =  G0id  ~  constant  *  E  F  •  S  (C-l) 
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The  term  G  represents  any  of  the  quantities  R,  RU,  RV,  RW,  E,  and  the  term  F  represents 
the  fluxes  correponding  to  these  quantities  (see  Section  3.0).  The  quantity  F  •  S  is  referred  to 
as  the  flux  accumulator. 

Suppose  G  =  R.  The  flux  accumulator  for  R  is  the  variable  DR.  The  constant  in  Eq. 
(C-l)  involves  At  and  volume,  i.e,  DT  and  VOL,  and  is  called  DTSVOL.  So,  if  we  write  Eq. 
(C-l)  above  for  G  =  R,  we  obtain 

R(JX,JY,JZ)  =  R(JX, JY, JZ)  -  DTSVOL(JX)  *  DR(JX,JY,JZ)  (C-2) 

which  is  precisely  one  of  the  equations  in  subroutine  UPDATE.  This  equation  updates  the 
old  value  of  R  by  the  net  mass  which  has  entered  the  volume  during  the  time  interval  At, 
hence,  the  subroutine  name. 

Now,  Eq.  (C-l)  represents  either  part  of  Eq.  (7),  Ref.  1;  i.e.,  the  predictor  and  the 
corrector  steps  have  the  same  form.  Equation  (C-2)  above  is  valid  for  either  step  also.  Which 
step  is  actually  being  computed  at  some  instant  in  the  program  depends  on  a  complicated 
manipulation  of  indices  and  counters  which  takes  place  in  several  subroutines.  These 
manipulations  are  beyond  the  scope  of  the  present  discussion,  other  than  to  observe  that 
they  occur. 

6.  DTCAL 


This  subroutine  determines  the  time  step,  At,  to  be  used  in  the  predictor-corrector 
calculations.  It  would  be  desirable  for  this  time  step  to  be  as  large  as  possible  so  that  a 
solution  could  be  obtained  quickly  and,  thus,  economically.  However,  the  computational 
scheme  is  unstable  for  time  steps  greater  than  some  critical  value.  By  unstable  we  mean  the 
solution  diverges  as  time  increases  rather  than  converges  as  we  desire.  The  critical  value  is 
determined  by  the  Courant-Friedrichs-Lewy  stability  criterion*,  which  says  that  the 


movimnm  n  llrv\troV\!o  +Jma  r  +  id  +Viq  mitiitvmm  Af 

muAimum  ctiiw  vvuuiv  lxjliiv  oi-wp  u  mv  iiiiiiiiiium  ui 


throughout  the  mesh,  where  V  is  the  volume  of  an  element,  S  and  q  are  area  and  velocity 
vectors,  respectively,  and  c  is  the  local  speed  of  sound.  This  criterion  essentially  says  that  the 


*Anderson,  J.  D.  Modern  Compressible  Flow.  McGraw-Hill,  New  York,  1982,  pp.  312-313. 
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time  step  can  be  no  greater  than  the  least  time  it  takes  for  an  acoustic  signal  to  travel  from 
one  volume  center  to  the  next,  which  is  certainly  physically  reasonable. 


DTCAL  proceeds  by  computing  at  a  given  time 


|q  ♦  S|  +  c | S | 


for  every  volume  in  the  mesh  and  searching  for  the  minimum  value.  It  is  interesting  and 
instructive  to  look  at  a  few  details  of  how  this  quantity  (called  DT)  is  actually  calculated. 

Careful  study  of  DTCAL  shows  that,  in  the  notation  of  this  program  but  leaving  off  the 
subscripts,  DT  is  determined  by  the  expression 


DT  = 


VOL  *  R 


( 


RQDA  |  +  I  RQDRQ  *  ADA  * 


GS2 


E/P  -  RGM1 


) 


'/2 


where 


RQDA  =  gq  «A 
RQDRQ  =  gq  •  gq  =  g2 q2 
ADA  =  A  •  A  =  A2 

and  GS2/[(E/P)  -  RGM1]  =  1/M2,  if  the  definitions  of  GS2,  RGM1,  E,  and  P  are  used. 
This  expression  for  DT  can  then  be  written  as 


DT  = 


G  *  VOL 


6|q  •  A|  +  (G2q2A2/M2) 


Vt 


But  M  =  q/c,  so  this  reduces  to 


DT 


VOL 

q  •  A  |  4-  cA 


which  is  the  same  as  the  equation  for  At. 

7.  FLUX 

The  subroutines  FLUXX,  FLUXY,  and  FLUXZ  (collectively  referred  to  as  FLUX) 
compute  the  fluxes  of  the  various  flow  variables  across  each  pair  of  surfaces  of  a  volume 
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element.  Thus,  FLUXX  computes  the  flux  across  the  two  A1  surfaces,  FLUXY  computes 
the  A2  fluxes,  and  FLUXZ  determines  the  transport  across  the  two  A3  surfaces.  These 
routines  also  increment  the  flux  accumulators,  which  symbolically  are  F  •  S,  with 
appropriate  care  for  the  signs.  The  three  routines  are  identical  in  form  and  differ  only  in  the 
subscripts  and  variables  used. 

The  computations  in  FLUX  are  straightforward.  Central  to  these  calculations  is  the 
definition  of  the  flux  of  something: 


flux  =  (  (  SOmethin8  ^  q  .  dA 
\  unit  volume  / 
s  x  ' 

which,  for  planar  surfaces  with  constant  flow  properties,  reduces  to 


The  “something”  here  is  the  dimensionless  mass,  (X,  Y,  Z)  momentum,  or  energy,  and  is 
simply  R,  RU,  RV,  RW,  and  E.  The  area  used  in  the  dot  product  depends  on  the  subroutine 
being  used.  Thus, 

q  •  A  =  QDA  =  (RU*AX1  +  RV*AY1  +  RW*AZ1)/R  in  FLUXX, 

=  (RU*AX2  +  RV*AY2  +  RW*AZ2)/R  in  FLUXY,  and 
=  (RU*AX3  +  RV*AY3  +  RW*AZ3)/R  in  FLUXZ 

The  equations  in  FLUX  then  follow  from  combinations  of  the  “something  per  unit  volume” 
and  the  ODA  terms. 

What  is  not  so  straightforward  in  FLUX  (and  in  some  of  the  other  subroutines  as  well)  is 
the  way  in  which  the  indices  (or  subscripts)  and  counters  are  manipulated.  For  instance,  in 
FLUXX  we  find  the  following  indices  and  DO  LOOP  parameters  being  used:  KODE,  JJ, 
JX1,  JX2,  NEX,  NNY,  NNZ,  JS,  JX,  JY,  JZ,  and  JT.  FLUXY  and  FLUXZ  have  similar 
quantities.  Actually,  this  rather  complicated  situation  is  the  result  of  economizing  the 
computations.  The  solution  process  involves  performing  a  set  of  calculations,  namely  those 
in  Eq.  (11),  Ref.  1,  for  every  volume  element  in  the  computational  mesh  and  at  every  time 
step.  Within  Eq.  (11)  we  see  that  there  is  a  need  for  repeated  determinations  of  F  •  S, 
although  the  subscripts  and  superscripts  appearing  in  Eq.  (11)  suggest  that  different  values 
of  F  and  S  are  to  be  used  in  the  various  terms.  This  is  indeed  true,  although  with  careful 
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study  of  the  overall  computational  process  we  find  that  certain  calculations  are  repeated 
exactly.  For  example,  at  a  given  time  step  the  quantity  F  •  S£  for  one  volume  element  is 
exactly  the  same  (except  for  a  sign)  as  F“  •  for  a  neighboring  volume  element. 
(Reference  1  and  Section  4.0  may  help  clarify  this.)  In  any  case,  the  important  idea  is  that 
certain  types  of  calculations,  namely  F  •  S,  are  done  very  many  times  in  this  program  and 
some  are  repeated  exactly.  The  design  of  the  FLUX  subroutines  is  such  that  the  equations 
for  mass,  momentum,  and  energy  flux  can  be  used  for  every  volume,  every  time  step,  and 
each  form  of  F  •  S  in  Eq.  (11).  The  versatility  of  these  subroutines  comes  from  the 
manipulation  of  the  indices  and  counters. 

A  number  of  subroutines  besides  FLUX  are  involved  in  these  manipulations;  however, 
much  of  this  activity  does  take  place  in  FLUX,  and  a  feel  for  what  goes  on  can  be  obtained 
by  studying  some  of  the  details  in  the  FLUX  routines,  which  we  shall  now  do  (on  a  very 
limited  scale). 

FLUX  evaluates  the  quantity  F  •  S.  The  particular  values  of  F  and  S  used  in  any  given 
calculation  depend  on  which  volume  is  being  considered,  which  face  of  the  the  volume  is 
being  considered,  the  time  step,  and  which  step  of  the  predictor-corrector  sequence  is  being 
executed.  Suppose  we  fix  the  volume,  face,  and  time  step  in  question.  A  representative 
situation  is  sketched  below. 
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The  volume  in  question  is  the  middle  volume,  the  face  of  interest  is  the  shaded  one  labelled 
SJ.  With  the  time,  volume,  and  face  fixed,  we  only  need  to  worry  about  which  step  in  the 
predictor-corrector  cycle  we  are  in.  The  two  steps  in  the  predictor-corrector  cycle  are 
represented  by  the  two  equations  in  Eq.  (11),  Ref.  1.  Looking  at  these  equations  we  see  that 
we  use  F  °  S£  in  one  step  and  F+  <>  Sf  in  the  other.  (Don’t  worry  about  the  bar  over  F  + 
right  now.)  So,  in  one  step  of  the  predictor-corrector,  we  use  the  value  of  F  “behind”  the 
face  in  question;  in  the  other  step  we  use  the  value  of  F  “ahead”.  This  can  be  accomplished 
by  using  the  same  equation  for  F  •  S  and  simply  changing  the  subscript  of  F  in  the  k 
direction.  (Remember,  the  subscripts  of  a  flow  quantity  refer  to  the  volume  in  which  it 
occurs.)  If  the  k  direction,  in  this  example,  is  in  the  X  direction,  then  we  increment  JX  to  go 
from  one  side  of  the  face  to  the  other.  Such  a  procedure  must  be  done  for  all  six  faces  of  this 
volume.  The  face  labelled  S^in  the  sketch  would  use  F_  and  F,  and  so  on. 

Besides  determining  F  °  S  at  each  face,  FLUX  sums  or  accumulates  the  net  flux  into  each 
volume  by  incrementing  the  flux  accumulators  DR,  DRU,  DRY,  DRW,  and  DE.  However, 
I  care  must  be  taken  with  the  signs.  For  instance,  in  the  sketch,  the  product  F  •  S£  might 

represent  the  flux  out  of  the  middle  volume  or  the  flux  into  the  right-hand  volume.  For  one 
volume  this  number  is  positive,  for  the  other  it  is  negative.  The  subscripts  and  signs  must  be 
carefully  matched  to  properly  sum  or  accumulate  the  fluxes. 

The  foregoing  discussion  has  established  the  reasons  for  the  multitude  of  indices  and 
index  operations  in  the  program.  The  following  brief  summary  shows  how  some  of  the 
subscripts  provide  a  means  to  switch  from  one  side  of  a  surface  to  another  and  to  switch 
from  one  volume  to  another. 

Each  of  FLUXX,  FLUXY,  and  FLUXZ  has  a  dimensioned  variable  KODE  (JJ)  whose 
16  values  are  specified  in  a  DATA  statement  in  each  subroutine.  The  value  of  JJ  is 
established  in  STEP  and  is  transferred  through  the  labelled  COMMON  CONST.  JJ  tells 
'  FLUX  which  phase  of  the  predictor-corrector  sequence  to  compute  by  picking  a  value  of 

KODE.  KODE,  in  turn,  sets  JX1,  JY1,  or  JZ1  depending  on  which  FLUX  subroutine  we  are 
using.  JX1  (or  JY1  or  JZ1),  in  turn,  sets  the  value  of  JT.  JT  sets  the  value  of  JS,  and  JS  tells 
FLUX  whether  to  use  the  flow  value  on  one  side  or  the  other  of  a  surface  in  evaluating 
F  •  S  and  whether  to  give  this  quantity  a  plus  or  minus  sign  when  it  is  added  to  the  flux 
accumulator  (schematically  represented  here  by  D). 

Figure  C-4  presents  the  values  of  and  relations  between  the  various  indices  as  the 
calculations  in  FLUX  proceed.  Also  depicted  are  representative  volume  elements  and 
reference  surfaces  to  illustrate  the  predictor-corrector  selection  and  flux  summation. 
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JX1  =  1  +  JT  =  2,  JS  =  JX,  D(JS)  =  D(JS)  -  F(JS-l)  +  F(JS); 
q  upstream  of  A1 

JX1  =  2  -*■  JT  =  3,  JS  =  JX-1,  D( JS)  =  D( JS )  -  F(JS)  +  F(JS+1); 
q  downstream  of  A1 

a.  FLUXX 


JYl  =  1  JT  =  2,  JS  =  JY,  D( JS )  =  D( JS)  -  F(JS-l)  +  F(JS); 
q  below  A2 

JYl  =  2  +  JT  «  3,  JS  ”  JY-1,  D(JS)  =  D(JS)  -  F(JS+1); 
q  above  A2 

b.  FLUXY 


JZl  =  1  JT  =  2,  JS  =  JZ,  D( JS)  =  D( JS)  -  F(JS-l)  +  F(JS); 
q  counterclockwise  of  A3 

JZl  =  2  -*•  JT  =  3,  JS  =  JZ-1,  D( JS )  =  D( JS )  -  F(JS)  +  F(JS+1); 
q  clockwise  of  A3 

c.  FLUXZ 

Figure  C-4.  FLUX  calculations  and  representations. 
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